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Abstract. Detailed calculations of the transport coefficients of a recently 
introduced particle-based model for fluid dynamics with a non-ideal equation of 
state are presented. Excluded volume interactions are modeled by means of biased 
stochastic multiparticle collisions which depend on the local velocities and densities. 
Momentum and energy are exactly conserved locally. A general scheme to derive 
transport coefhcients for such biased, velocity dependent collision rules is developed. 
Analytic expressions for the self-diffusion coefficient and the shear viscosity are 
obtained, and very good agreement is found with numerical results at small and 
large mean free paths. The viscosity turns out to be proportional to the square root 
of temperature, as in a real gas. In addition, the theoretical framework is applied 
to a two-component version of the model, and expressions for the viscosity and the 
difference in diffusion of the two species are given. 



PACS numbers: 47.11. +j, 05.40. +j, 02.70.Ns 
1. Introduction 

The interplay between hydrodynamic interactions and thermal fluctuations is crucial 
for a wide range of phenomena in biophysics and soft-matter physics. Several particle- 
based simulation methods such as dissipative particle dynamics [H [2] , smooth particle 
hydrodynamics [3] , and direct simulation Monte Carlo [H [5] have been developed for 
the efficient modeling of these phenomena with the motivation to coarse-grain out 
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irrelevant atomistic details while correctly incorporating the essential physics. One 
particular method was introduced by Malevanets and Kapral in 1999 [6l [7j and is often 
called stochastic rotation dynamics (SRD) [H [9l [TOl [TTl [12] or multi-particle collision 
dynamics (MPC) |l3l [lU |15l |16]. The method is based on so-called fluid particles 
with continuous positions and velocities which follow a simple, artiflcial dynamics and 
undergo efficient multi-particle collisions. The original SRD algorithm models a fluid 
with an ideal gas equation of state. Recently, the collision rules of the original algorithm 
have been generalized to model excluded-volume effects, allowing for a more realistic 
modeling of dense gases and immiscible binary fluids [TTl [H] • These new multi-particle 
collision rules depend on local densities and velocities. The new models can be thought 
of as a coarse-grained multi-particle collision generalization of a hard sphere fluid since, 
just as for hard spheres, there is no internal energy. Their static properties such as the 
equation of state and the entropy density have been derived recently [TTl HH] , and were 
shown to be in excellent agreement with numerical results for the pressure, speed of 
sound, density fluctuations and the phase diagram. However, the dynamic properties 
are not yet fully understood. Some preliminary results about the transport coefficients 
were published previously [TTl [THl [191 [20] , but without derivation and without results 
for the viscosity at small mean free path. In this paper I present the details of how 
to derive expressions for the self-diffusion coefficient and the shear viscosity at small 
and large mean free path. For simpler, unbiased collison rules, this was done before 
[HI [TOl [m [211 [22] ; the added difficulty here is that the collision probabihties depend on 
the particle velocities which leads to strong correlations. The derived formulas will be 
compared to numerical simulations in various limits. 

The paper is structured as follows: section 2 introduces the simulation model, and in 
section 3 the self-diffusion coefficient is derived. In section 4 an equilibrium method 
to determine the kinetic part of the viscosity is used: Green-Kubo relations of the 
kinetic stress correlations are evaluated. Section 5 applies a non-equlibrium method 
to calculate the coUisional contribution to the viscosity by evaluating the momentum 
transfer in shear flow. This non-equilibrium method is a more rigorous version of the 
earlier approach introduced in ref. [8] and now, in principle, can be used to obtain 
the shear dependence of the coUisional viscosity. In section 6 the developed scheme is 
applied to calculate transport coefficients for a binary model and a model with a large 
collision rate. 
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2. Model 



As in the original SRD algorithm, the solvent is modeled by a large number N of point- 
like particles of mass Trip which move in continuous space with a continuous distribution 
of velocities. The particle mass is set equal to one throughout this paper. The system 
is coarse-grained into {L/aY cells of a d-dimensional cubic lattice of linear dimension 
L and lattice constant a. The algorithm consists of individual streaming and collision 
steps. In the free-steaming step, the coordinates, rj(t), of the solvent particles at time 
t are updated according to Vi{t + t) = Vi{t) + T\'i{t), where Vi(t) is the velocity of 
particle i at time t and r is the value of the discretized time step. In order to define 
the collision, we introduce a second grid with sides of length 2a which (in d = 2) groups 
four adjacent cells into one "supercell". 

As proposed in Refs. [8[[9], a random shift of the particle coordinates before the collision 
step is required to ensure Galilean invariance. All particles are therefore shifted by the 
same random vector with components in the interval [—a, a] before the collision step 
(Because of the supercell structure, this is a larger interval than in the conventional 
SRD algorithm). Particles are then shifted back by the same amount after the collision. 
To initiate a collision, pairs of cells in every supercell are randomly selected. As shown 
in fig. 1, three different choices are possible: a) horizontal (with <ti = x), b) vertical 
{(T2 = y), and c) diagonal collisions (with cr^ = {x + y)/V^ and (7^ = (x — y)/\/2). 
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Note that diagonal colhsions are essential to equilibrate the kinetic energies in the x— 
and y— directions. 

In every cell, we define the mean particle velocity, 

Un = -TT- / Vj, (1) 

i=l 

where the sum runs over all particles, M^, in the cell with index n. The projection of 
the difference of the mean velocities of the selected cell-pairs on crj, Au = aj ■ (ui — U2), 
is then used to determine the probability of collision. If Au < 0, no collision will be 
performed. For positive Am, a collision will occur with an acceptance probability which 
depends on Am and the number of particles in the two cells. Mi and M2. This rule 
mimics a hard-sphere collision on a coarse-grained level: For Am > clouds of particles 
collide and exchange momenta. The following acceptance probability was introduced 
in ref. [T7| 

Pa(Mi,M2, Am) = 6'(Am) tanh(A) with A = v4AmMiM2, (2) 

where 6 is the unit step function and A is a parameter which allows us to tune the 
equation of state. The hyperbolic tangent function was chosen in (|2]) in order to obtain 
a probability which varies smoothly between and 1. This paper treats the two most 
relevant limits for the acceptance probability: In the first one, which was proven to 
be thermodynamically consistent, one keeps the collision parameter A small such that 
the tank can be replaced by its argument. In the other limit, pa is simply the theta 
function alone which is equivalent to let A go to infinity. In principle, any other 
situation with A < 00 can be handled by expanding the tank; all additional terms still 
contain solvable integrals. 

Once it is decided to perform a collision, an explicit form for the momentum transfer 
between the two cells is needed. The collision should conserve the total momentum 
and kinetic energy of the cell-pairs participating in the collision, and in analogy to 
the hard-sphere liquid, the collision should primarily transfer the component of the 
momentum which is parallel to the connecting vector cTj. This component will be called 
the parallel or longitudinal momentum. There are many different choices which fulfill 
these conditions. Since the original goal was to obtain a large speed of sound, a collision 
rule was selected which leads to the maximum transfer of the parallel component of the 
momentum and does not change the transverse momentum. The rule is quite simple; 
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it exchanges the parallel component of the mean velocities of the two cells, which is 
equivalent to a "reflection" of the relative velocities, 

v!{t + r)-J = -{v!{t)-J), (3) 

where m" is the parallel component of the mean velocity of the particles of both cells. 
The perpendicular component remains unchanged, 

vHt + r)=v^{t). (4) 

More specifically, the rule for horizonal collisions in cri direction is simply 

with Vi^x denotes Vi^xii + T)- The y-components of the particle velocities stay the same, 
Vi,y = Vi^y. For the upward diagonal collision along the 0-3 direction one has 

Vi,y = Ux + Uy- Vi^x (6) 

It is easy to verify that these rules conserve momentum and energy in the cell pairs. 

Because of x — y symmetry, the probabilities for choosing cell pairs in the x— and y— 
directions (with unit vectors cri and (T2, see fig. 1) are equal, and will be denoted 
by w. The probability for choosing diagonal pairs {cr^ and 0-4, see fig. 1) is given 
by Wrf = 1 — 2ty. w and Wd must be chosen such that the hydrodynamic equations 
are isotropic and do not depend on the orientation of the underlying grid. This was 
done by considering the temporal evolution of the lowest moments of the velocity 
distribution function, or alternatively, by evaluating the isotropy of the kinetic part 
of the viscous stress tensor [181 ISO]- Both approaches concluded that w = 1/4 and 
Wd = 1/2, and previous simulations show that both the speed of sound and the shear 
viscosity are isotropic for this choice. Note, however, that this does not guarantee that 
all properties of the model are isotropic. This becomes apparent at high densities or 
high collision frequency, 1 /r ^ 1 , where the cubic anisotropy of the grid shows leading 
to inhomogeneous states with cubic symmetry [17J . 
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3. Diffusion coefficient 



The self-diffusion constant D is given by a sum over the velocity-autocorrelation 
function (eq. (102) in [9]), 

oo ' 

D = rY, {vMMnT)), (7) 

n=0 

where the prime on the sum indicates that the n = term has the relative weight 1/2. 
The particle mass is set to one throughout the paper. Assuming molecular chaos, one 
has 

{v^{0) v^{nT)} = keT g where g = — (8) 

and the sum ([7]) turns into a geometric series which is easily summed up to give 
knTr 1 + q 

D = — — . 9 

2 1-g 

g is the velocity correlation one time step apart, which completely defines the diffusion 
coefficient in the molecular chaos approximation. There is four different contributions 
to as a result of possible horizontal, vertical, and the two diagonal collisions: 

g = w{gH + gv) + Y^3u + 9D) (10) 

The vertical, gv, and the horizontal contribution, gn, occur with probability w = 1/4, 
whereas the upward diagonal part, gu and the downward diagonal part, go, enter 
with probability Wd = 1/2. The additional factor 1/2 in front of the second bracket 
accounts for the fact that the two diagonal collisions are always chosen together and 
on average only 1/2 of the particles in a 2a x 2a-super cell undergo an upward diagonal 
collision while the other half is subject to a downward diagonal collision. In contrast, 
in a horizontal (or vertical) collision operation all particles in the super-cell can in 
principle be involved in the collision since the super cell is divided in two double cells 
with horizontal (or vertical) collisions in each one. The contribution due to vertical 
collisions is easy to calculate, since in this case, the value of Vi^x is unchanged in a 
collision, eq. (jl]), and Vi^xi^) = '^j,x(0). Following eq. we find 

= k^ = ' 
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Consider now the horizontal contribution, gn- We assume that there are n particles in 
the left subcell, and m particles in the right subcell of the cell pair where the collision 
happens. It is also assumed, that the tagged particle i is contained in the left subcell. 
For the moment the particle numbers n and m are kept fixed. According to eq. ([2]), 
for small A, a collision occurs with the acceptance probability, 

Pa = Anm 9{Au) Au] (12) 

and the value of Vi^^ changes in the following way Vi^xi^) = — t'j^^(O). Here, is 
the x-component of the mean velocity of both subcells: 



^ n+m 

Ux = — : Vj^cc (13) 

n + m ^ 
i=i 

For horizontal directions. Am contains only the x-components of the particle velocities: 

Am = - V Vj^x V Vj^^ . (14) 

n ^-^ m ^-^ 

j=l j=n+l 



If a collision is not accepted, which happens with probability 1 — pa, the value of Vi^^ 
remains unchanged: Vi^xi^) = ^i,x(0). From these considerations it follows, 

9H = -j^{pAVi,xi2Ux - Vi^x) + {l-pA)ViJ . (15) 

The average is performed over the velocity distribution (at time t = 0) of all involved 
particles and over the particle numbers in the two subcells forming a collision cell. 
Note that there is a non-trivial coupling between the velocity-dependent acceptance 
probability pa and the updated velocity t>j .^(r). Ignoring this coupling, which amounts 
to the approximation (p^ .r(0)fj_^(r)) ^ {pa) {vj^x{0)vj^x{T)) leads to an errorous result 
which is a factor of about two smaller than the correct one, even for large particle 
numbers n, m. 

Using (vf^x) — ^bT, eq. (fT5l) simplifies to 

2 2 

5-// = 1 + T^{PAVi,x{,Ux - Vi^x)) = 1 + -r^\Knm)n,m (16) 
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The average {■■■)n,m denotes an average over the particle number distribution in the 
two subcells (which are Poisson-distributed and uncorreletad in an ideal gas). Kn^n 
contains the average of the velocity distribution, which is Maxwellian if molecular 
chaos is assumed; and it has to be split up in two parts, 

Knm — 2 {^nm ^nm) (1'^) 

because there are two possibilities occuring with probability 1/2, respectively : (i) the 
tagged particle i can be in the left subcell as one of the n particles there, or (ii) can be 
in the right subcell with a total of m particles. Hence, one has 

Km^^nmJl^ with 

/+0O n+m 
dv^+^ e (Ax.) I^u v,,,{u, - v^,,) n h{vi,.) ■ (18) 
-~ i=i 

fo{vi,x) is the Maxwell-Boltzmann distribution. 



'l,X 



Uvix) = Texp , with T = , . (19) 

Assuming the tagged particle to be in the left subcell in the definition of J^^ means 
also that n > in the expression for J^^. For the same reason, m must be larger than 
zero in expressions for J^^. 

Because of the ^-function the direct evaluation of the {n + m)-dimensional integral 
becomes very tedious for large n, m. Therefore, we express the ^-function by means of 
the (5-f unction, 

poo 

e{Au) = / 5{Au-c) dc (20) 
Jo 

which allows us to use the integral-representation of the 5-function, 

exp{ik{Au - c)) ^ (21) 
-oo 27r 

Even though this transformation leads to two more integrations it enables a 
straightforward calculation of J^^ for arbitrary n and m. It involves only simple 
integrals of type exp{—x^) with n = 0, 1, 2, 3. The integral to be solved is: 

/•oo poo 11 r+oo 

jL = r"+- dc ^ dv?t"^ 



nm 



27r 
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n ^ n+m \ ^ 

j=l j=n+l / 

n ^ n+m \ / ^ n+m 

j=l j=n+l / \ j=l 



2 I _|_ 2 \ 



f22) 



First, the argument of the exponential has to be transformed into a sum of squares in 
order to decouple the velocities and to enable integrating over them. This is achieved 
by the transformations for the velocities in the left subcell: 



ksTik . 

Wi = Vix for I = l...n 

n 

and for the ones in the right subcell, 

ksTik 



(23) 



Wi = Vi,x + 



for i = n + l...n + m . 



m 



(24) 



Expressing eq. fl22l) in the new variables Wi gives: 

""^ dk '•+°° 



/' = r 



n-\-m 



dc 







27r 



^^n+m 



/ ., kBT-fk^\ 1 r 2 2 2 l\ 

X exp \-tkc 1 exp ( Fl + ^2 + ••• + ^n+mj 1 



Wi + W2 + .-+Wn Wn+l+Wn+2 + ■■■+U!n+r 



n 



m 



+ ksTik'j 



Wi + 



kBTik\ f wi + W2 + ■■■ + w, 



n 



n+m 



n + m 



Wi 



ksTik 



n 



(25) 



with 7 = {l/n) + (1/m). The integration over the Wi is straightforward and yields 

Tn+m 



r = V 



. , dk 
dc I — 







X exp —ike — 



^ 27r 
keT^k"' 



S 



2 

-+7 



n 



1 



1 



(26) 

n + m 

The remaining integral over h and k is done easily using similiar transformations. The 
final result is: 



3/2 



n^ 



1 1 
- + — 

n m 



1 - 



1 



n + m 



(27) 



which is always negative. 
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The quantity J^^ is given by an expression almost identical to eq. ([27]) , one merely has 
to replace vi^x by Vn+i,x at the two positions where vi^x appears isolated. This changes 
the index of the tagged particle from 1 to n + 1, meaning that this particle is not in 
the left sub cell but in the right one. One finds = leading to an expression for 
Jnm which is symmetric in n and m: 



Anm 



nm ' nml 



AjkBT) 



3/2 



2nm 



+ n + m — 1 



(2J 



Using eqs. (11611181) gu is calculated. In an ideal gas, Knm has to be averaged over 
the poisson-distributed particle numbers in the subcells, n and m. In the current 
model with a non-ideal equation of state, particle number fluctuations are supressed 
compared to an ideal gas. Therefore, we will neglect these fluctuations completely 
for the moment, and replace n and m by the mean number of particles in a subcell 
M = (n) = (m) . This leads to the final result 



9h 



1-2A 



^^3/2 



(29) 



3.2. Diagonal collisions 

Here, the upward diagonal collision (along as in fig. [1]) is considered first. In this 
case, the velocity difference. Am, depends on both the x- and the y-component of the 
particle velocities. 



Am = CTg ■ (Ui - U2) = {Ui^x - U2,x + Ui^y - U2,y) 



(30) 



with the normal vector of the diagonal collision direction, 0-3 = (1, l)/-\/2, leading to 
a more involved calculation of gu. Furthermore, both components of the velocities 
change in a collision, but only the change of the x-component is needed: 



Vi,xir) = Ux + Uy - Vi 



(31) 



This change occurs with probability p^, see eq. ( IT2l) . but with Am given in eq. ( l30l) . 
It follows 



1 



9u 



{pa Vi,x{Ux + Uy- Vi^y) + (1 - Pa)vIx) ■ 



(32) 
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The average must be done over the velocity distribution (at time t — 0) oi all involved 
particles and over the particle number in the two subcells forming a collision or double- 
cell. Because of (vf ^) — ksT we can write 

9U = 1 + 7 rp i-^nm ~ (33) 

with the abbrevations 

Rnm = {pa Vi,x{Ux + Uy - Vi^y))T (34) 
Lnm^{PAvl^)T (35) 

The subscript T means that only the thermal average is taken; the particle numbers n 
and m of the subcells are kept fixed in this average. 

We consider again the two equivalent situations that the tagged particle is in the left 
subcell with a total of n cells denoted by the superindex /, or in the right cell, described 
by the superindex r. Both possibilities occur with equal probability 1/2 and we have 

Lnm = - {Lnm + ^nm) (36) 

Other abbreviations are introduced as Rnm = ^nmFnm, and Lnm = Anmlnm- This 
time the y-componcnt of the velocities shows up in the calculations, leading to a 
2n + 2m-dimcnsional integral for F^„^. Wc again adopt the trick to express the 9- 
function in the acceptance probability p^ by means of an integral over the 5-function 
and use the integral representation of d{x), which yields the following integral: 

roo poo jr, r+oc r+oo 

Km = r"+- / dc - dv-+- / d^J- 

Jo J-oc J -oo J -oo 

{/ 1 [l " 1 1 \ 

I 7f (^^.^ + '"^'V^ ~ m ^ ^^''^ ^ '"''"^ ~ ^ I 

\ L j=l j=n+l J / 

X exp {-^^ (^M + - + ^n+m,x + + - + ^n+m,J } 
n ^ n+m \ 

j=l J=n+1 / 

n+m \ 



X — 

V2 \n 



X 



n + m 
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This integral is solved straightforwardly in a very similar fashion to the one in eq. 
The final result is 



{ksT) 



3/2 



V2^ 

A similar calculation leads to 



1 

nm 



2n2 



(3J 



7 + 



2n2 



(39) 



Again, we observe the symmetry that the quantities with upper index r follow from 
the ones with index I by simply interchanging n and m in the obtained expressions: 
= -F^n, I^m = ^Inn- Inserting these results into eq. ( |33ll and neglecting particle 
number fluctuations gives 



9u 



I- A 



M3/2 



(40) 



The evaluation of the downward diagonal contribution was done in a similar way, and 
for symmetry reasons one finds go = gu- According to eqs. fl29|) . fHOl) and ([9]), ffTOj) 
the diffusion coefficient follows in the limit of small collision parameter A as 



D = kBTT\^ 



keT 2 



(41) 



which is in good agreement with simulation data, see fig. 2. It is interesting to 
note, that the contributions to D from the vertical and horizontal collisions together, 
{gn + fi'y)/4, are exactly equal to the ones from the diagonal ones, {gu + gn)/'^, which 
is a sign of the correct isotropic choice of the probabilities w and Wd- 



3. 3. Density fluctuations 

Density fluctuations were completely ignored in the derivation of eq. 0411) . In principle, 
the correct density fluctuations which follow from the equation of state could be 
incorporated, but here we will only estimate the effect of these fluctuations. This 
can be done by assuming that the particle numbers n and m in the two subcells are 
uncorrelated and behave like in an ideal gas. Instead of replacing n and m by the mean 
value M in the expressions for gu and gu we average over a Poisson-distribution. In 
this case the probabihty to flnd n particles in the left subcell is equal to 



Pceu,n = exp(-M) M"'/n\ , 



(42) 
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where M is the average particle number per cell, M = (n). On the other hand, any 
one of these n particles can be our tagged particle of index i. Hence, the probability 
that a given particle is in a subcell with a total of n particles is equal to Pn ~ npceii,n 
(see for instance Ref. [H]). Normalizing this probability gives 

1 e-*^M" , , 

Care must be taken for the case when at least one of n or m are zero. The formal 
expression might look divergent, but one has to keep in mind how the algorithm is 
actually working: The acceptance probability pa ~ nm is zero in this case, hence no 
collision is executed. This means that the velocities do not change, fi,x(T) = fi,a:(0), 
and the actual value of gn and gu is one. This can be incorporated by starting the 
summations over n and mat?T, = l,m=l instead of zero in eq. fj44|) . 

There is another subtle point about this average: For quantities with index Z, such as 
^nmy '^^ have to use Pn to average over the particle number in the left subcell which 
contains the tagged particle, but Pceii to average over the particle number m in the 
right subcell. This is because no specific particle is adressed in the right subcell; all 
what is needed is the probability to find a given particle number m in this cell. These 
arguments lead to the following result. 



, //crT / / nm \m n 

Qh = 1- 2A\ ( J \ hn + m-1 

^ \ 2tt \\l n + m l2n 2m 



+ PmPcell 



KRi \ ^ \ ^ / nm 



r\ / , / , \ PnPcelLm\ I 

ivr ^ — ^ ^ — ^ I y n + m 

n=l m=l ^ ' 

n 1 

1 — (n + m — 1] 

2m 2^ 



ml 

1 — (n + m — 1) 

2n 2^ ^ 



nm 



^ . ... . , (44) 

n + m 

Similar expressions hold for the diagonal terms, gu and go- These averages over 
the Poisson-distributed particle numbers cannot be done analytically. Therefore, a 
numerical evaluation of the sums was performed and inserted in expressions ( iQfTOj) 
to obtain the diffusion coefficient. This result provides an upper limit for the effect of 
fluctuations. A check shows that it only differs by a few percent from the result without 
fluctuations. As seen in figure 2 of ref. [IT], the simulation results show excellent 
agreement with the formula without particle number fiuctuations. The effective 
excluded volume interactions of the current non-ideal model lead to a supression 
of density fiuctuations. Hence, it seems plausible that neglecting particle number 
fiuctuations entirely is a better approximation for this model than assuming strong 
ideal-gas like density fiuctuations. 
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The multi-particle collision algorithm consists of two steps - streaming and collision. 
Each of these steps redistribute momentum. At large mean free path A compared to 
the cell size a, this transport mainly occurs in the streaming step where momentum 
is directly carried away by the moving particles. At small mean free path, most 
momentum is transfered in the collision step. The momentum transport during 
streaming is characterized by the so-called kinetic viscosity, whereas the transfer during 
collisions gives rise to the coUisional contribution to the viscosity. The kinetic viscosity 
is given by the Green-Kubo relation, [2T] . 



lykin = rj^ Y^ {<yxy,k{f^)<yxy,k{nr)) , (45) 

^ n=0 

The prime on the sum indicates that the first = term has a prefactor 1/2. a^y^k is 
the off-diagonal element of the kinetic stress tensor, 

N 



In order to evaluate eq. (I46ll the correlation of the stress tensor for only one time step 
apart is needed, since under the assumption of molecular chaos, the sum in eq. (H5ll is 
again a geometric series. Hence, the following quantity is required, 

N 

b = i^iA^)^iA'^)H^i^)Hyi^)) ■ (47) 

Analogous to eq. ( ITOl) we have 

b = w{bv + bH) + ^{bu + bn) (48) 

where by describes the contribution from the vertical collisions, Bh the horizontal, and 
bu and b^ the ones from the diagonal collisions. 



4.1. Horizontal contribution 

During a horizontal collision, only Vj^x can change. Thus, using {viaVjp) = SijSap ksT 
we have 

N 

bn = Y {[PAVi,xVi,y[2ux - Vj^x]vj^y + (1 - PA)vi,xVi,yVj,xVjJ) (49) 
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N{kBTf + 2kBT(pAj2 




(50) 



where the time argument was omitted since all velocities are now at the same time 
t = 0. The sum over all particles can be rewritten as a sum over all double cells and 
over the particles in each of them since both Ux and pa depend only on the properties 
of the particles in the considered double cell. This leads to the following simplification: 



Nh = N/{2M) is the number of double cells for horizontal collisions. The sum runs 
over the n particles in the left subcell of a double cell and over the m particles in the 
right subcell. The average has to be taken over n, m and over the velocity distribution. 
Inserting the definition of p^, eq. (fT^ . yields 



As before, the subindex T denotes the thermal average at fixed particle numbers per 
cell, and the subindex n, m refers to the average over the particle numbers. 

The remarkable feature of eq. (1521) is that the calculation of bn is reduced to the same 
integrals which appeared in the previous section in the calculation of the diffusion 
coefficient, namely the quantities ^nm defined in eqs. ( IT61IT71) . This also 

means as in the case of the diffusion coefficient, there again is a non-trivial coupling 
between the velocity-dependent acceptance probability and the elements of the stress 
tensor. As before, a mean-field like treatment consisting of multiplying the averaged 
value of Pa with the averaged stress-correlations leads to a large error - a factor of two. 
This error persists even in the case of high particle numbers where one naivly would 
assume this decoupling to be correct. According to eq. fll7|52l) we have 




(51) 




(52) 




(53) 



Using the explicit form of the K' s, eqs. ( fTTfTSl) . gives 




(54) 
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Ignoring particle number fluctuations, and replacing n and m by their mean value M 
leads to the final expression: 



bn = N^ksTf I 1 - 2A\ ^M^/' , \ (55) 



TT 



which looks identical to gn, eq. (129|) . except the prefactor. 

Since the stress tensor is symmetric with respect to interchanging x and y, and 
the collision rules are constructed to be also x-y-symmetric, one has for the vertical 
contribution: by = bn- 



4-2. Diagonal contributions 

The upward diagonal contribution can be written as 

N 

bu = ^ {[PAVi,xVi,yVj^^Vj^y + (1 - PA)Vi,xVi,yVj,xVj,y]) (56) 
N 

= N{kBTf + ^ {pAVi,xVi,y[vj,xVj,y - Vj^^Vj^y]) (57) 

The tilde denotes velocities at time r as opposed to time zero. The sum over the particle 
index j is rewritten as a sum over all available double cells for diagonal collisions and 
another sum over the particles in that cell. The reason again is that both the velocities 
at time r and the collision probability pa depend only on the considered double cell. 
One finds, 

I f N \ n+m \ 

bo = N{kBTf + Nr,( I ^ Vi^^Vi^y j ^ PA[vj,xVj,y - Vj^xVjJ ) (58) 

where = N/ (4M) is the number of double cells for diagonal collisions. (This 
number is smaller than the one for horizontal collisions due to the specific algorithm 
chosen) Inserting the collision rules for upward diagonal collisions, eq. [H]), yields 

bo = NiksTf + 

I f N \ n+m \ 

'^y '^j,x '^j,y ] (59) 

AN I f ^ \ I 

= N{kBTf + v,^xv,,^y^ {nm ^ 9iAu) Au (60) 

x{ux + Uy)[ux + Uy- Vj^x - fiJ)r>„,^ (61) 
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Ua and Am are properties of the collision cell and can be put in front of the sum. The 
summation gives zero, since _ = Q. Thus, 



For the downward diagonal collisions a similar calculation results in bo = hjj for 
symmetry reasons. 

Final result for kinetic viscosity 

Combining the results for by = bn, eq. ( l55l) . and bo = bu, eq. (|62l) . gives 



With r = b/ {N{kBTy') the summation of the geometric series leads to the kinetic 
viscosity 



which is exactly equal to the diffusion coefficient, eq. (l^T]) and also scales with v^fc^T. 
Hence, in the limit of large mean free path the Schmidt number, Sc = ly/D, of this 
model goes to one. Sc can be modified by adding additional regular SRD rotations on 
the subcell level which do not change the equation of state. 

5. Collisional viscosity 

In contrast to the kinetic viscosity, the collisional contribution is evaluated in non- 
equilibrium — in shear fiow with shear rate 7. Only the limit of infinitesimally small 
shear rates will be considered, that means the collisional viscosity will be calculated to 
zeroth order in 7. In principle, it is possible to perform the calculation in higher order 
which would tell us about eventual shear- or shear thinning similar to ref. [22l [23] . 
however, in lowest order a number of convenient approximations can be made. For 
example, we can neglect the deviation of the particle velocity distributions from a 
Gaussian shape since this distortion is caused by the very small shear rate and gives 
rise to a higher order term in the viscosity. Due to the shear, the velocity distribution 



bu = N{kBTf 



(62) 




(63) 




(64) 
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depends on the height yi of particle i and the following factorized form of the N-particle 
distribution function can be assumed, 

N 

fN,x ~ JJ foivi,x - iVi) (65) 

i=l 
N 

fN,y ~ n^O^^*'?') ^66) 

i=l 

where /o is the Maxwell-Boltzmann distribution, eq. f|T9l) . This factorization 
is equivalent to assuming Molecular chaos which was shown to be a very good 
approximation for calculations of the coUisional part of transport coefficients [21]. A 
given collision cell (which consists of two sub cells) is divided by a line y = h, and 
average transfer of x-momentum across this line is calculated. Vertical collisions do not 
change the x-component of the velocities at all, and hence do not have to be considered. 
However, both diagonal and the horizontal collisions transfer x-momentum across this 
line and will be evaluated in the next subsections. 



5.1. Horizontal collisions 

A dividing line is assumed to go across a pair of horizontally aligned subcells at fixed 
height y = h with < h < a, see fig. El The left subcell contains A^^ particles which 
are split into two populations below and above the dividing line with particle numbers 
n, m, respectively, such that A^^ = n + m. The line divides the particle population in 
the right subcell into p and q particles, with total particle number Nb = p + q. For 
now, the particle numbers n, m, p, and q are kept fixed. I consider the change of total 
x-momentum of the n + p particles below the dividing line, which equals the amount 
of x-momentum transferred over this line to the particles denoted as filled circles in 
fig. [21 This momentum transfer Ap^ for fixed height h and fixed particle numbers is 
averaged over the particle positions and velocities resulting in 

(Ap.)Lp,= ^^p^5f Ai;,,.^ ^ . (67) 
The upper index, H, stands for "horizontal" and the averages are defined as: 

POO 

{■■■)v= / -/^..^r^^^" (68) 



oo 



{■■■)y = TL; dyr \^ 1.^^, / -dyr' (69) 



h 
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Later, additional averages over the particle numbers and the height h (which is 
equivalent to an average over the random shifts) will be performed. The dashed line in 
fig. [2] divides the two cells into four cell fractions; the restricted position average, (l69l) . 
takes into account that particles can be at any height y within their corresponding 
cell fractions with equal probability. No averaging over the lateral positions is needed 
since the flow speed depends on y but not on x. For small collision parameter A the 
acceptance probability is given by eq. f|T2|) and is proportional to 9{Au). Au is equal 
to the difference of the mean velocities in the two subcells projected onto the cri = x 
direction, 

^ n+m ^ n+m+p+q 

Au = Ua - Ub = — ^ Vi,x - ■ C^'O) 

^ i=l ^ i=n+m+l 

The particle number in the left subcell is Na = n + m, and A''^ = p + q is the 
particle number in the right subcell. In order to simplify the integration over the 
6'-function, the exponential representation from eqs. fl20|21l) is utilized again. Avi^x 
is the change of x-component of a particle and follows from the collision rule (JSj) as 
Avi^x = '^{ux — Vi^x)- Here, Ux is the center of mass velocity of the two subcells together, 
Ux = {vi^x + ■■■ + Vn+m+p+q,x)/{NA + Nb)- Eq. ([67]) bccomes 

nmpq 

2ANaNb (^(^J^^ dcj^ ge^fc(A-c)^^ ^(^^ _ \^ (71^ 

The velocity average is performed over the approximated non-equilibrium velocity 
distribution /at^^., defined in eq. dHSD which still depends on the height yi of the particles. 
This suggests the following variable transformation: 

Vi,x = Vi^x - iVi ■ (72) 

In addition, the auxiliary variable c is changed into c by the transformation 

c = c - 7(5 

which moves the |/j-dependence from the exponent into the integral boundaries. One 
obtains, 

{^Px)nmpq ~ 

I t-CO POO 7, poo 

2ANaNbT'''^"'+p+i ( dd — exp[ik{Au - d)] 

\J-^Q J-oo J-oo 
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n+p 

[{n + p)u^-J2i^j,x + iR])Y (74) 



with 



n+m+p+q n+p 

In order to decouple the velocity integrals, the argument of the exponential is 
transformed into a sum of squares by changing the velocity variables in the left subcell, 

Wi = Vi^x — for i^l...NA (76) 

and the ones in the right subcell, 

Wi = Vi,x + -^^ iori = Na + I...Na + Nb. (77) 

It is straightforward to perform the resulting Gaussian integrals over the new variables 
Wi, however, it is sufficient to only evaluate the term linear in the shear rate 7. There is 
also a shear independent term, but it will cancel out later in the final average over the 
dividing line. This is to be expected since there should be no net momentum transfer 
in equilibrium where 7 = 0. The next higher order term must be cubic in the shear 
rate, since the viscosity (which is proportional to the momentum transfer divided by 
the shear rate) should not depend on the sign of the shear rate. The relevant linear 
term has a simple structure, 

(Ap.>Lp. = 0(1) + 0(f ) + 

2ANaNb 7 keT (^^ dc ^exp {^-ikc - ^jML^ ^fc , ^ (78) 

^ = Q(&-Tr)+R^ (79) 



Na Nb, 

One can show that replacing the lower shear dependent boundary in the c-integral by 
zero will only cause an error of order 7^ in the momentum transfer which is negligible. 
With this modification the integrals over k and c have simple solutions, and one arrives 
at 



(Ap.)L,, = ANANBiyl^ {Q)y + 0(1) + 0(f) (80) 
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Figure 2. A double cell consisting of two subcells for liorizonal collisions along the 
(Ti direction defined in fig. [1] ua is the horizontal component of the mean velocity 
of the n + m = Na particles in the left subcell. The horizontal component of the 
p + q ~ Nb particles in the right subcell is given hy ub- 



The position average (...)y defined in eq. f l69l) over particle heights give (?/i)y = h/2 
for particles below the dividing line with i = l,...n. Above the dividing line, where 
i = Na + 1, ■■■,Na + p, one finds (?/j)y = (a + h)/2. These results help finding the 
following averages, 

a f m q 

a {n + p){m + q) 
2' 



{Q)y 

{R)y 



Y 



Na + Nb 
2{pm + nq) + nm + pq 



pq_ 



nm 



NaNb 

The numbers in the cell fractions, n and m, as well as p and q are fiuctuating even 
at fixed Na and Nb and are approximately binomially distributed. This leads to the 
following rules when performing the binomial average {■■■)b over the particle numbers, 

=NaW, I^p)^ = NbW 

{m)b =Na{1-w), {q)b = Nb{1-w) 

{nm)b = {Na - l)w{l - w) , {pq)b = {Nb - l)w{l - w) , (82) 

where w = h/a is the relative height of the dividing line. Using eqs. (I80ti82p and 
averaging over the height of the dividing line, {...)h = (l/c^) Jq ...dh, the average 
momentum transfer at fixed A^i and Nb is obtained. 



Nb 



A'ja IkBT 



{Na - Nb 



\21 



NaN^ 



B 



+ 0(7 



^3) 
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5.2. Diagonal collisions 

Both diagonal collisions along 0-3 and cr^ in fig. 1 will also transfer x-momentum over 
a horizonal diagonal line. Only the contribution from the upward diagonal, 0-3, will 
be discussed here, since the contribution from the down-wards diagonal gives the same 
result for symmetry reasons. It is neccessary to distinguish between the two cases 
depicted in figs. [3] and HI where (A) the dividing hne goes through the lower cell and, 
(B) where the line goes through the upper cell. Only both contributions together will 
give the correct result with the same structure than eq. fl83l) . 

First, case A is considered. The change of x-momentum of the particles described by 
full circles in fig. [3] is given as 

{ApXmp=(^(pAj2Av,,}j ^ (84) 

where the upper index. A, denotes upward diagonal collisions for situation A. The 
average over the particle height and the velocity distribution differs from 

the previous definition eqs. fl68|69p 

POO 

(••>= / •••/^,./^,.rf<r^^rf<r^^ (85) 



00 



pa /•2a 



the velocity average now involves the y-component of the particle velocities as well. 

The acceptance probability pa is again given by eq. (fT2l) and proportional to 6{A.u) 
which however has now a different form, 

^ n+m ^ n+m+p 

^.,^.,,-u, = -j:^^^-^ -^. (87) 

1=1 i=n+m+l 

The particle number in the bottom left cell is Na = n + m, and A^^ = p, is the particle 
number in the top right cell. Using the exponential representation (120112 ip for the 6- 
function, and the upward colhsion rule ([6]) to redefine Avi^x = Ux + Uy — Vi^x — Vi^y, one 
finds 

{Apx)^^pg = ANaNb 



X 



' POO POO IT \ 



V Y 
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This expression is evaluated using similar substitutions and approximations than for 
horizontal collisions. After performing all averages one obtains the momentum transfer 
for fixed particle numbers Na_, Nb (in linear order in the shear rate 7): 



T2"V 2^ 



AT, 



(iV^ _ 1) ( 1 _ __£ ) + uNaNb 



^9) 



A similar calculation for case B where the dividing line cuts the top right cell, see fig. 
m results in 



A'ja ksT 
T2"V 2^ 



{Nb - 1) 1 



Na 
Nb 



UNaNb 



(90) 



By averaging both results, one finds the momentum transfer for the upward (U) 
diagonal collision for A^i > 1 and Nb > 1: 

1 
2 



{Na 



knT 



6 



27r7 



^NaNb + 



^NaNb 



(91) 



This expression is symmetric in Na and Nb-, as it should be, and has the same general 
structure than eq. ( !83l) . Finally, this result is multiplied by two to account for the 
contribution from the downward diagonal and added together with the result for 
horizontal collisions using the weight factors w and Wd for horizontal and diagonal 
collisions, respectively. The total average momentum transfer is: 

{{^Px))na,Nb = 

{Na-Nb^ 



Aja / kBT 
^2~V 2^ 



{28wd + 8w)NaNb + {wd + 2w)- 



(92) 



NaNb 

Now, one can derive the collisional part to the kinematic viscosity, I'coii, which is 
proportional to the momentum transfered per time and length 



2a T 



7 t^coU P 



(93) 



p = M/ a^ is the density of our two-dimensional fluid. Replacing the particle numbers 
Na and Nb by the average particle number per cell, M and setting w = 1/4 and 
Wd = 1/2, one gets 



(94) 
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for not too small mean particle number per cell, M > 2. Figs. O and E] compare the 
theoretical result for the total viscosity, u = VcoU + J-'kin from eqs. (16411941) with numerical 
results obtained by evaluating Green-Kubo relations for stress autocorrelation functions 
(for numerical details see refs. [IZIEO]). One sees that the agreement at small and large 
mean free paths is very good. The small deviations at intermediate mean free path are 
within the error bars of the data. However, there is another possibility; in regular SRD 
it was found, [21], that at A <^ a corrections due to the violation of molecular chaos 
become relevant in i/^j^, but not in Ucou- These effects are large for very small A, but 
in this limit they are irelevant for the total viscosity because it is dominated by Vcoii- 
Hence, if at all, one would see these deviations only at intermediate mean free paths. 
It is possible that there are similar deviations for the current model which could show 
at intermediate A. 



There is a simple argument for the structure of expression (]94l) : Collisions "smear 
out" momentum over a distance of order cell size a in one time step r. The kinematic 
viscosity is the coefficient of momentum diffusion which in analogy to a random walk 
is defined by the "hopping " distance a per time step as -D ~ o^/t. However, the 
collisions occur with some collision rate which is proportional to the thermal average 
of the acceptance probability given by 

{pa) = A {e{Au) AuNaNb) ~ A y/k^ M^/^ (95) 

Therefore, one would expect 

i^coii ~ -{pa) ~ - a ^/k^M'/' (96) 
r T 

which is exactly what was found. A similar argument can be made for the kinetic part 
of the viscosity. The analogy to a random walk predicts Ukm ~ A^/r since momentum 
is now only "hopping" a distance of order mean free path A = Ty^ksT. In addition, it 
is clear that at — there is almost no collisions, the momentum is travelling huge 
distances and the kinetic viscosity becomes very large. This can be described by an 
additional factor 1/ {pa) and one obtains. 



1 A 



. , ^~^v^M-3/^ (97) 

{Pa) t {pa) a 

which, apart from a small correction term, agrees with what was derived rigorously, 
eq. f lMl) . Using this conceptual insight, one can easily predict the general form of 
these expressions for arbitrary acceptance probabilities without actually going through 
tedious derivations. 
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6. Results for other models 



6.1. Large collision rate, A ^ oo 



The calculation of the previous sections are easily modified to treat the limit of large 
collision parameter A oo, where the acceptance probability is simply pA = 6{Au). 
One finds, 



2M - 1 + e 



for the kinetic viscosity if ideal gas-like fiuctuations are assumed. The exponential 
terms come from averaging over the particle number fiuctuations assuming a Poisson- 
distribution. This assumption is correct for an ideal gas, but not for the current model 
with a non-ideal equation of state, which has smaller fiuctuations. Furthermore, it 
was shown that the density fiuctuations in the limit of large A are even smaller than 
predicted by the equation of state [19j. Given the small size of the exponential terms 
~ g-2A/ ^Qj. typical particle numbers M between three and twenty, the fiuctuation 
corrections, i.e. the exponential terms, can be safely neglected. 

By modifying the non-equilibrium calculation scheme of section 5, one obtains the 
average momenum transfer in shear fiow for Na > 1, Nb > 1 as 

{{APx))NA,Nn = 

[{Uwd + 4w)NANB + {wd + 2w){NA + NB-2)] (99) 



24{Na + Nb 

Ignoring particle number fiuctuations by setting A^^ = A''^ = M, substituting w = 1/4, 
Wd = 1/2, the coUisional viscosity follows from eq. fl93|) . 



coll 



4M V M 



(100) 



Fig. [7] shows excellent agreement of simulation data for M = 3 from ref. [r9J with the 
theoretical results for the total viscosity u = Ukm + i^coii from eqs. ( l98l llOOp . One sees 
that the viscosity is basically independent of the particle number M at small mean free 
path. 



6.2. Transport coefficients for a binary mixture 



In [18] a multi-component version of the current collision model was introduced. 
Consider a binary system with two types of particles, A and B. In order to obtain 
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phase separation, a repulsion between different kind of the particles was implemented; 
but there is no repulsion among particles of the same kind. This is done in the following 
way: Suppose a double cell is selected for a possible collision. A-particles in cell 1 can 
undergo collisions with B-particles from cell 2. For symmetry reasons, B-particles 
from cell 1 are also checked for possible collision with A-particles in cell 2. The rules 
and probabilities for these collisions are exactly the same as in the one-component 
situation. This means, most of the results derived above can be used with only minor 
modifications. 

Ma = (Na) is the average number of A-particles in a cell of size a, Mb = {Nb) is 
the average of B-particles in such a cell. The total particle density in a cell is now 
given by p = {Ma + MB)/a^. While using eqs. f|T6|) to fl27j) for diffusion calculations 
one has to distinguish whether the tagged particle is of the A or B type. For the 
self- diffusion of A-particles one has to set n = Na and m = Nb in eq. (ITHi) for K^^, 
but n = Nb and m = A^^ in -ftT^^ which, because of the symmetry = K[^^, leads 
to Knm = K^NaNb- "^^^ Corresponding results for the auxiliary variables qh, Qu and 
g are not symmetric with respect to the particle numbers Na and Nb-, and a short 
calculation gives g^ = ^ — 20"^, g^ = 1 — 0^, and g^ = 1 — (p^ with 



and 7 = I/Ma + ^/Mb where the actual numbers A^^, Nb were already replaced by 
their average values. Ma and Mb, respectively. Similar considerations hold for the 
diffusion of B particles and one obtains the diffusion coefficient for A and B particles, 
respectively. 



Both particles diffuse the same way only if Ma = Mb- In this case the diffusion 
coefficient D = Da = Db agrees with the one for the single- component model, eq. 
( 14T|) . If there is less A than B particles, the diffusion of A particles is slower even if the 
masses of both particles are the same as assumed here. This is physically plausible, 
since a single A-particle in a sea of B particles will be backscattered very often whereas 
the many B-particles keep going and are less affected by the collisions. Fig. [8] plots 




(101) 




(102) 
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Figure 3. A double cell consisting of two subcells for collision along the upward 
diagonal direction cr^ as defined in fig. [1] ua is the component of the mean velocity 
projected on cr^, of the n + m = Na particles in the left subcell. The projected mean 
velocity of the p — Nb particles in the upper right subcell is given by ub- A line 
divides the lower cell at height h, < h < a. 



the ratio of the two diffusion coefficients Da/ Db as a function of the averaged relative 
particle number difference, Ap = {Ma — Mb) /{Ma + Mb) for fixed total number 
M = Ma + Mb- Ap was only varied in a range were both particle numbers Ma and 
Mb are always larger or equal to one because density fluctuations become significant 
at smaller particle numbers but are neglected in the theory. For a total number of 
M = 5 and Ma = 1, Mb = 4 one sees that B-particles diffuse about twice as fast than 
A-particles. 

Finally, the kinetic and collisional viscosities are determined by applying concepts from 
the previous sections, and one finds 



bin 
kin 



bin 
coll 




ksTMAMB 
3r V 2tt{Ma + Mb) 



AMaMb + 



{Ma 



Mb)' 



AMaMb 



(103) 



7. Conclusion 



In this paper, I have presented a detailed, systematic derivation of the transport 
coefficients of a multi-particle collision model for fluid flow with a non-ideal equation 
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Figure 4. A double cell consisting of two subcells for collision along the upward 
diagonal direction as defined in fig. [1] ua is the component of the mean velocity 
projected on 0-3 of the p — Na particles in the lower left subcell. The projected mean 
velocity of the n + m — Nb particles in the upper right subcell is given hy ub- A line 
divides the upper cell at height h, a < h < 2a. 
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Figure 5. Total viscosity as a function of time step r for small collision rate. 
The symbols are numerical results obtained by evaluating Green-Kubo relations for 
M = 3. The solid line is the theoretical result, v = Vcoii + Vkin from eqs. (|64|94p . 
Parameters: L = 64a, M = 5, fcsT = 1, and A = 1/60. 
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Figure 6. Same as fig. [5j total viscosity v = Vcoii + Vkin as a function of time step 
T, but for larger r, and compared to the theoretical result. 




time step 



Figure 7. Total viscosity v = Vcoii + Vkin as a function of time step r for large 
collision rate, A — > oo. The symbols are numerical results obtained by evaluating 
Green-Kubo relations for M = 3. The solid line is the theoretical result from eqs. 
(|98|100p . for M = 3. The dashed line is the theoretical result for AI — 10. Parameters: 
L = 64a, ksT = 1, and A oo. 



of state. Similar calculations have been published before, but to my knowledge, this is 
the first derivation which accounts for velocity-biased collision rules in MPC Analytic 
expressions for the self- diffusion coefficient, and the shear viscosity are obtained, and 
very good agreement is found with numerical results at small and large mean free paths. 
For the analysis of the collisional contribution to the viscosity, an earlier scheme was 
improved which can now be used to derive possible shear thinning or thickening. The 
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Figure 8. Ratio of tlie diffusion coefficient for A and B particles, Da/Db, as a 
function of averaged relative particle number Ap = {Ma — Mb)/M for fixed total 
number M = Ma + Mb- The dashed line is the theoretical result from eqs. (|102p . for 
M = 3. The solid line is for M — 5, and the dotted one is for M — 10. Parameters: 
kBT =1, 1/60, T 1, Ma > 1, Mb > 1. 

general formalism was also applied to a binary collision model with a miscibility gap. 
General arguments for the scaling of the transport coefficients with temperature and 
particle number as a function of the collision rule are given. This allows one to predict 
without tedious derivations how the transport coefficients change if the collision rule 
is modified. 

For the current collision rule, the viscosity turns out to be proportional to the square 
root of temperature as in a real gas. It is shown that at large mean free path, the 
viscosity and the diffusion coefficient are exactly equal, resulting in a Schmidt number 
of one in this limit. 
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